library(ggplot2)
library(list)

###load survey
load("JanCourier.RData")

#################Variable information
####All variables with a political figure's name indicate support for the figure (1=support, 0=not support)

####Contemporary represents the number of contemporary figures supported on a list, which is the combined 
#vector representing the number of figures supported on the list from both those respondents asked only the 
#number of figures they support (qv2_3 and qv2_4) and the number of figures they support and do not support 
#(qv2_1a and qv2_2a)

####History represents the number of historical figures supported on a list, which is the combined 
#vector representing the number of figures supported on the list from both those respondents asked only the 
#number of figures they support (qv1_3 and qv1_4) and the number of figures they support and do not support 
#(qv1_1a and qv1_2a)

###THistory and TContemporary are vectors representing whether or not an individual was in a control or 
#treatment group for the respective experiments (1=Treatment, 0=Control)


###################Analyses

##direct question regarding Putin support, Table 1
mean(na.omit(jan$Putin))

#########################################################
#### Analyses of relative support for contemporary figures
##########################################################

####frequency table, Table A1
table(jan$Contemporary,jan$TContemporary)

###examine experimental results for different types of counting
####Not reported, but described in fn 18
####list both support and not support, support answer
t.test(jan$qv2_1a,jan$qv2_2a)
###list just support question
t.test(jan$qv2_3,jan$qv2_4)
###test if different methods of counting yield different counts
t.test(jan$qv2_1a,jan$qv2_3)
t.test(jan$qv2_2a,jan$qv2_4)

##test for overall treatment effect, Table 1
t.test(jan$Contemporary[which(jan$TContemporary==0)],jan$Contemporary[which(jan$TContemporary==1)])

###dummy vector for nas
cont.na<-jan$Contemporary[is.na(jan$Contemporary)==F]
Treatment.na<-jan$TContemporary[is.na(jan$Contemporary)==F]

###Blair/Imai test for design effect, fn 25
ict.test(cont.na,Treatment.na,3)

#####Vectors for political figures directly supported
mean(na.omit(jan$Ziuganov))
mean(na.omit(jan$Zhirinovsky))
mean(na.omit(jan$Mironov))

########Vector of direct support for contemporary figures
jan$CS2<-jan$Ziuganov + jan$Zhirinovsky + jan$Mironov


###floor effects
######## support Putin and no others, fn 15
mean(na.omit(jan$Putin == 1 & jan$Zhirinovsky == 0 & jan$Ziuganov == 0 & jan$Mironov==0))

###at 1 for contemporary treatment, pg 11
mean(na.omit(jan$Contemporary[jan$TContemporary==1])==1)

###modify contemporary experiment so that individuals at one in treatment receive 0, rerun t-test, fn 22
C2<-jan$Contemporary
C2[which(jan$Contemporary==1 & jan$TContemporary==1)]<-0
t.test(C2[jan$TContemporary==0],C2[jan$TContemporary==1])

####mean number in control, list then direct; not reported, akin to fn 12
mean(na.omit(jan$Contemporary[which(is.na(jan$CS2)==F & jan$TContemporary==0)]))
mean(na.omit(jan$CS2[which(is.na(jan$Contemporary)==F & jan$TContemporary==0)]))

###regression of deflation, Table A2
summary(lm(Contemporary ~ TContemporary*CS2, data=jan))

###plot the relationship between direct responses and experimental responses, figure 2
##turn treatment into factor
Treatment<-as.factor(jan$TContemporary)
levels(Treatment)<-c("Control","Treatment")

jan$Condition<-Treatment
qplot(CS2,Contemporary, data=jan, xlab="Number of control-group figures supported in direct questions", 
      ylab="Number of figures supported in list", alpha=I(.0001)) +
  geom_smooth(aes(linetype=Condition), color="black", size=1) + theme_bw() + theme(legend.position=c(.2,.8)) +
  scale_linetype_manual(values=c("solid","dotted")) 

###Get predictions of SDB, Blair/Imai method Table 1
set.seed(123)
summary(dim.r <- ictreg(Contemporary ~ 1, data = jan, treat = "TContemporary", J=3, method = "lm"))
pdim.r<-predict.ictreg(dim.r)

df<- glm(Putin ~ 1, data= jan, family=binomial)

predict(dim.r, direct.glm=df)

#################################
####Repeat with historical data
##################################

###examine differences based on wording
####Not reported, but described in fn 18
###support and not support question, support answer
t.test(jan$qv1_2a,jan$qv1_1a)

###just support question
t.test(jan$qv1_4,jan$qv1_3)

###test for difference in counts based on counting technique
t.test(jan$qv1_1a,jan$qv1_3)
t.test(jan$qv1_2a,jan$qv1_4)

####frequency table, Table A1
table(jan$History,jan$THistory)

###test for overall treatment effect, Table 1
t.test(jan$History[which(jan$THistory==0)],jan$History[which(jan$THistory==1)])

###remove na from vector
history.na<-jan$History[-1080]
Treatment.na<-jan$THistory[-1080]

###evidence of design effect, using Blair and Imai method, fn 25
ict.test(history.na,Treatment.na,3)


####support for historical figures, pg 9
mean(na.omit(jan$Stalin))
mean(na.omit(jan$Brezhnev))
mean(na.omit(jan$Yeltsin))

###vector of direct support for all historical figures
jan$HS2<- jan$Stalin + jan$Brezhnev + jan$Yeltsin 

####################Floor effects
########Support for Putin and no others, fn 15
mean(na.omit(jan$Putin == 1 & jan$Stalin == 0 & jan$Brezhnev == 0 & jan$Yeltsin==0))

###at 1 for history treatment, pg 10
mean(na.omit(jan$History[jan$THistory==1])==1)

###modify historical experiment so that individuals at one in treatment receive 0, rerun t-test, fn 22
H2<-jan$History
H2[which(jan$History==1 & jan$THistory==1)]<-0
t.test(H2[jan$THistory==0],H2[jan$THistory==1])

### difference between individuals in control in number supported, nas removed; experiment first; not reported, akin to fn 26
mean(na.omit(jan$History[which(is.na(jan$HS2)==F & jan$THistory==0)]))
mean(na.omit(jan$HS2[which(is.na(jan$History)==F & jan$THistory==0)]))

####deflation regression, table A2
summary(lm(History ~ THistory*HS2, data=jan))

####Analyze relationship betweeen number of control figures supported and experimental results, Figure 2
Treatment<-as.factor(jan$THistory)
levels(Treatment)<-c("Control","Treatment")

jan$Condition<-Treatment
qplot(HS2,History, data=jan, xlab="Number of control-group figures supported in direct questions", 
ylab="Number of figures supported in list", alpha=I(.0001)) +
  geom_smooth(aes(linetype=Condition), color="black", size=1) + theme_bw() + theme(legend.position=c(.2,.8)) +
  scale_linetype_manual(values=c("solid","dotted")) 

###Get predictions of SDB, Blair/Imai method, Table 1,
set.seed(123)
summary(dim.r <- ictreg(History ~ 1, data = jan, treat = "THistory", J=3, method = "lm"))
pdim.r<-predict.ictreg(dim.r)

df<- glm(Putin ~ 1, data= jan, family=binomial)

predict(dim.r, direct.glm=df)

############################################
################Double experiment analysis, Table 1
#########################################
A<-cbind(jan$History[jan$THistory==1],jan$Contemporary[jan$THistory==1]) ##Create matrix for treatment
B<-cbind(jan$History[jan$THistory==0],jan$Contemporary[jan$THistory==0]) ##create matix for control

A<-A[!is.na(A[,1]), ] ##remove NAS lw
A<-A[!is.na(A[,2]), ] ##remove NAs lw (no NAs in control)

###Average estimated support
(mean(A[,1])- mean(B[,1])+mean(B[,2])- mean(A[,2]))/2
####### Glynn formula for CI
(mean(A[,1])- mean(B[,1])+mean(B[,2])- mean(A[,2]))/2+qnorm(.975)*sqrt((1/(2*(nrow(A) + nrow(B))))*(var(A[,1]) + var(A[,2]) - 2*cov(A[,2],A[,1]))+(1/(2*(nrow(A) + nrow(B))))*(var(B[,1]) + var(B[,2]) - 2*cov(B[,2],B[,1])))
(mean(A[,1])- mean(B[,1])+mean(B[,2])- mean(A[,2]))/2-qnorm(.975)*sqrt((1/(2*(nrow(A) + nrow(B))))*(var(A[,1]) + var(A[,2]) - 2*cov(A[,2],A[,1]))+(1/(2*(nrow(A) + nrow(B))))*(var(B[,1]) + var(B[,2]) - 2*cov(B[,2],B[,1])))

######################################
#### Analyze correlation between figures, 
########################################

####Analyze all correlations, footnote 14
cor(na.omit(jan[,c(9,12:14,17:21,22:25)]))

